/*
 *  lib-smv: Matrix - Vector Library targeted to Structural Mechanics Problems
 *  Copyright (C) 2006-2008  St.Brcic, Lj.Zugic
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
/** @file Cholesky.cpp
 *  @brief Implementation of class Cholesky - LU decomposition of SPD matrices
 *
 *        Standard decomposition     A = L * LT 
 *        Generalized decomposition  A = L * D * LT
 *
 *  DIRECT SOLUTION OF LINEAR EQUATIONS A * x = b or A * X = B
 *  
 *  LU decomposition of symmetric positive definite (SPD) matrices used for:
 *  (1) Solution of linear equations
 *  (2) Inverse of a SPD matrix
 *  (3) Determinant of a SPD matrix
 * 
 */ 

#include "Cholesky.h"
namespace smv_lib
{

//! empty (default) constructor
Cholesky::Cholesky(void)
{
  isMatrix_A    = false;
  isNameGiven   = false;
  isVector_b    = false;
  isMatrix_B    = false;
  isStandard    = false;
  isGeneralized = false;
  isDecomposed  = false;
  isInverse     = false;
  isSolved      = false;
  isMore        = false;
}

//! constructor with a given matrix A
//! and a choice of decomposition method (default is standard: А=L*LT)
Cholesky::Cholesky(const Matrix& AA, DECOMPOSE dec) : A(AA), Decompose(dec)
{
  N = A.getRows();
  assert(N == A.getCols());
  isMatrix_A  = true;
  Name = A.getName();
  if (Name == "MATRIX")
    isNameGiven = false;
  else
    isNameGiven = true;
  isVector_b  = false;
  isMatrix_B  = false;
  //--- check if A is symmetric and positive definite
  // is_SPD = checkInput(); 
  is_SPD = true; 
  setMethod(Decompose);
  isDecomposed  = false;
  isInverse     = false;
  isSolved      = false;
  isMore        = false;
}

//! constructor with a given matrix A and a vector b - solution of equations A*x=b
Cholesky::Cholesky(const Matrix& AA, const Vector& bb) : A(AA), b(bb)
{
  //-- matrix A
  N = A.getRows();
  assert(N == A.getCols());  // matrix A must be square
  isMatrix_A = true;
  Name = A.getName();
  if (Name == "MATRIX")
    isNameGiven = false;
  else
    isNameGiven = true;
  //--- right-hand side vector b
  assert (N == b.getLen());
  isVector_b = true;
  isMatrix_B = false;
  //--- solution vector x
  Vector* tmp = new Vector(N);
  assert(tmp != NULL);
  x = (*tmp);
  delete tmp;
  //--- check if A is symmetric and positive definite
  is_SPD = checkInput();
  setMethod(chol_STANDARD);
  isDecomposed  = false;
  isInverse     = false;
  isSolved      = false;
  isMore        = false;
}

//! constructor with a given matrix A and a matrix B - solution of equations A*X=B
Cholesky::Cholesky(const Matrix& AA, const Matrix& BB) : A(AA), B(BB)
{
  //--- matrix A (coefficient matrix)
  N = A.getRows();
  assert(N == A.getCols());   // matrix A must be square
  isMatrix_A = true;
  Name = A.getName();
  if (Name == "MATRIX")
    isNameGiven = false;
  else
    isNameGiven = true;
  //--- matrix B (right-hand side)
  assert (N == B.getRows());
  M = B.getCols();
  isMatrix_B = true;
  isVector_b = false;
  //--- matrix X (solution of equations A * X = B)
  Matrix* tmp = new Matrix(N,M);
  assert(tmp != NULL);
  X = (*tmp);
  delete tmp;
  //--- check if A is symmetric and positive definite
  is_SPD = checkInput();
  setMethod(chol_STANDARD);
  isDecomposed  = false;
  isInverse     = false;
  isSolved      = false;
  isMore        = true;
}

//! destructor
Cholesky::~Cholesky(void)
{
  if (Decompose == chol_STANDARD)
    delete LT1;
  if (Decompose == chol_GENERAL)
    delete LT2;
  if (isInverse == true) {
    delete L_inv;
    delete A_inv;
  }
}

//! set Matrix A
void Cholesky::setMatrix_A(const Matrix& AA, DECOMPOSE dec)
{
  //--- matrix A
  A = AA;
  N = A.getRows();
  assert(N == A.getCols());   // matrix A must be square
  isMatrix_A  = true;
  Name = A.getName();
  if (Name == "MATRIX")
    isNameGiven = false;
  else
    isNameGiven = true;
  //--- check if A is symmetric and positive definite
  is_SPD = checkInput();
  setMethod(dec);
  isDecomposed  = false;
  isInverse     = false;
  isSolved      = false;
  isMore        = false;
}

//! set the method of decomposition: standard or generalized
void Cholesky::setMethod(DECOMPOSE dec)
{
  Decompose = dec;
  if (Decompose == chol_STANDARD) {
    isStandard = true;
    isGeneralized = false;
  }
  if (Decompose == chol_GENERAL) {
    isStandard = false;
    isGeneralized = true;
  }
}

//! set the right-hand side vector b (Matrix A must be given first)
void Cholesky::setVector_b(const Vector& bb)
{
  if (isMatrix_A == false) {
    cerr << "\nGive the Matrix A first!";
    return;
  }
  b = bb;
  assert(N == b.getLen());
  isVector_b = true;
  //--- solution vector x
  Vector* tmp = new Vector(N);
  assert(tmp != NULL);
  x = (*tmp);
  delete tmp;
  setMethod(chol_STANDARD);
  isInverse     = false;
  isSolved      = false;
  isMore        = false;
}

//! set the right-hand side matrix B (Matrix A must be given first)
void Cholesky::setMatrix_B(const Matrix& BB)
{
  if (isMatrix_A == false) {
    cerr << "\nGive the Matrix A first!";
    return;
  }
  //--- matrix B (right-hand side)
  B = BB;
  assert (N == B.getRows());
  M = B.getCols();
  isMatrix_B = true;
  //--- matrix X (solution of equations A * X = B)
  Matrix* tmp = new Matrix(N,M);
  assert(tmp != NULL);
  X = (*tmp);
  delete tmp;
  setMethod(chol_STANDARD);
  isInverse     = false;
  isSolved      = false;
  isMore        = true;
}

//! check in the input matrix is symmetric and positive definite
bool Cholesky::checkInput()
{
  bool test = false;
  //--- check if A is symmetric
  isSymmetric = checkSymmetric();  // matrix MUST be symmetric to enter Cholesky class
  if (isSymmetric == false)
    cerr << "\nCholesky::checkInput():\n Matrix is NOT Symmetric!\n";
  //--- check if A is positive definite
  isPosDefinite = checkPosDefinite();
  if (isPosDefinite == false)
    cerr << "\nCholesky::checkInput():\n Matrix is NOT Positive Definite!\n";
  //--- symmetric and positive definite
  if (isSymmetric == true && isPosDefinite == true)
    test = true;
  //--- final result
  return test;
}

//! check if the matrix A is symmetric
bool Cholesky::checkSymmetric()
{
  if (isMatrix_A == false) {
    cerr << "\nCholesky::checkSymmetric():\n Matirx is not given!\n";
    return false;
  }
  double A_ij, A_ji;
  for (int i = 0; i < N; i++) {
    for (int j = (i+1); j < N; j++) {
      A_ij = this->A.getElement(i,j);
      A_ji = this->A.getElement(j,i);
      if (A_ij != A_ji)     // if at least one pair of symmetric elements are not equal
        return false;       // matrix is NOT symmetric
    }
  }
  return true;  // matrix is symmetrix
}

//! check if the matrix A is positive definite
bool Cholesky::checkPosDefinite()
{
  bool test = true;
  bool test1 = true, test2 = true;
  double a_ii, zero = 0.0;
  //---- 1st check: positive diagonal elements
  for (int ii=0; ii<N; ii++) {
    a_ii = this->A.getElement(ii,ii);
    if (a_ii <= zero)
      test1 = false;
  }
  //--- 2nd check: Sylvester's conditions:
  //    positive all minors related to diagonal elements
  //    check from Cholesky algorithm instead:
  int i,j,k;
  double sum, temp;
  //--- temporary matrices
  Matrix* A_tmp = new Matrix(A);
  assert(A_tmp != NULL);
  Matrix* L_tmp = new Matrix(N,N);
  assert(L_tmp != NULL);
  Vector* D_tmp  = new Vector(N);
  assert(D_tmp != NULL);
  //--- Colesky decomposition
  for (i=0; i<N; i++) {
    for (j=i; j<N; j++) {
      sum = A_tmp->getElement(i,j);
      for (k=0; k<i; k++)
        sum -= L_tmp->getElement(i,k) * L_tmp->getElement(j,k);
      if (i == j) {
        if (sum <= zero) {
          cerr << "\nCHOLDEC(): Matrix is NOT positive definite!" << endl;
          test2 = false;
        }
        D_tmp->setData(i,sqrt(sum));
      } // end if i == j
      else {
        temp = sum / D_tmp->getData(i);
        L_tmp->setElement(j,i, temp);
      }
    } // end for j=i
  } // end for i=0

  //--- combine checks: both must be true
  if (test1 == false || test2 == false)
    test = false;
  // release dynamic memory
  delete D_tmp;
  delete L_tmp;
  delete A_tmp;
  // final result
  return test;
}

//========================================== public working functions
//! Cholesky decomposition of symmetric and positive definite matrix
void Cholesky::Decomposition()
{
  //--- check if matrix A is a SPD matrix
  assert(is_SPD != false);
  //--- prepare Cholesky factors
  Matrix* tmp = new Matrix(N,N);
  assert(tmp != NULL);
  Vector* dd  = new Vector(N);
  assert(dd != NULL);
  //--- choose the method of decomposition
  switch (Decompose) {
      case chol_STANDARD:
         L1 = *tmp;
         D1 = *dd;
         cholDec();       // standard decomposition
         break;
      case chol_GENERAL:
         L2 = *tmp;
         D2 = *dd;
         cholDecGen();    // generalized decomposition
         break;
  }
  isDecomposed = true;
  //--- calculate also determinant of the matrix
  Det = Determinant();
  // release dynamic memory
  delete dd;
  delete tmp;
}

//! inverse of the matrix A
void Cholesky::Inverse()
{
  A_inv = new Matrix(N,N);
  assert(A_inv != NULL);
  //--- for inverse matrix use ONLY the standard decomposition (L*LT)
  if (isDecomposed == true && Decompose == chol_STANDARD) {
    *A_inv = cholInv();
    return;
  }
  //--- if previously general decomposition (L*D*LT) was called:
  //--- or, if NO decompostition was done, call the standard one
  else {
    //--- prepare room for matrices
    Matrix* tmp = new Matrix(N,N);
    assert(tmp != NULL);
    Vector* dd  = new Vector(N);
    assert(dd != NULL);
    L1 = *tmp;
    D1 = *dd;
    //--- now, call standard decomposition
    checkInput();
    cholDec();
    *A_inv = cholInv();
    // relese dynamic memory
    delete dd;
    delete tmp;
    return;
  }
}

//! determinant of the matrix A
const double Cholesky::Determinant()
{
  double di = 1.0, det = 1.0;
  int i;
  if (Decompose == chol_STANDARD) {
    for (i = 0; i < N; i++)
      di *= D1.getData(i);
    det = di*di;
  }
  else if (Decompose == chol_GENERAL) {
    for (i = 0; i < N; i++)
      di *= D2.getData(i);
    det = di;
  }
  return det;
}

//! solution of linear equations A*x = b or A*X = B
void Cholesky::Solution()
{
  if (isDecomposed == false) {
    cerr << "\nDecomposition must be done first!";
    return;
  }
  switch (isMore) {
      case false:
         x = cholEqs(b);    // single system of equations
         break;
      case true:
         X = cholEqs(B);    // more systems of equations
         break;
  }
  isSolved = true;
}

//! write results to a file (with predefined name "Cholesky.txt")
void Cholesky::Write()
{
  string FileName;
  if (isNameGiven == false)
    FileName = "Cholesky.txt";
  else
    FileName = Name + ".txt";
  this->WriteTo(FileName.c_str());
}

//! write results to a given file
void Cholesky::WriteTo(const char* file)
{
  //--- open the file
  out.open(file,ios::out);
  assert(out != NULL);
  out << "===================== CHOLESKY DECOMPOSITION ==========================" << endl << endl;
  out << "(1) INPUT MATRIX DATA: " << endl;
  out << "    - Name is:  " << Name;
  out << "    - Order N = " << N << endl << endl;
  out << A << endl;
  if (isDecomposed == true) {
    out << "(2) CHOLESKY FACTORS: " << endl;
    out << giveMethod() << endl;
    out << "Cholesky factor L" << endl;
    if (isStandard == true) {
      out << L1.getName() << endl;
      out << L1 << endl;
      out << "Cholesky factor LT" << endl;
      out << LT1->getName() << endl;
      out << (*LT1) << endl;
    }
    if (isGeneralized == true) {
      out << L2.getName() << endl;
      out << L2 << endl;
      out << "Cholesky factor D" << endl;
      out << D2 << endl;
      out << "Cholesky factor LT" << endl;
      out << LT2->getName() << endl;
      out << (*LT2) << endl;
    }
    out << "Determinant is " << Det << endl << endl;
    out << "(2a) CONTROL of decomposition: " << endl;
    out << "The product of Cholesky factors: " << endl;
    if (isStandard == true)
      out << " A = L * LT " << endl;
    if (isGeneralized == true)
      out << " A = L * D * LT " << endl;
    out << cholControl() << endl << endl;
  } // end if (isDecomposed == true)
  if (isInverse == true) {
    out << "(3) INVERSE MATRIX: " << endl << endl;
    out << (*A_inv) << endl;
    out << "Control: A * A_inv = " << endl;
    out << A * (*A_inv) << endl;
  }
  if (isSolved == true && isMore == false) {
    out << "(4) SOLUTION OF LINEAR EQUATIONS: " << endl << endl;
    out << "    Right-hand side = " << endl;
    out << b << endl;
    out << "    The solution is: " << endl;
    out << x << endl;
  }
  if (isSolved == true && isMore == true) {
    out << "(4) SOLUTION OF THE SYSTEM OF LINEAR EQUATIONS: " << endl << endl;
    for (int i=0; i<M; i++) {
      out << "   System number " << (i+1) << endl;
      out << "   Right-hand side = " << endl;
      out << B.getColData(i) << endl;
      out << "   The solution is: " << endl;
      out << X.getColData(i) << endl << endl;
    }
  }
  out.close();
}

//! give the method of decomposition as a string
const string Cholesky::giveMethod() const
{
  string info = "Method of decomposition is ";
  switch (Decompose) {
      case chol_STANDARD:
         info = info + "STANDARD: A = L * LT\n";
         break;
      case chol_GENERAL:
         info = info + "GENERALIZED: A = L * D * LT\n";
         break;
  }
  return (info);
}

//! get lower triangular Cholesky factor L
Matrix Cholesky::get_L()
{
  assert(isDecomposed == true);
  switch(Decompose) {
      case chol_STANDARD:
         return L1;
         break;
      case chol_GENERAL:
         return L2;
         break;
  }
  //--- if decomposition is not performed
  cerr << "\nDecomposition MUST be done first!\n";
  return A;
}

//! get transpose of the lower triangular Cholesky factor L, i.e upper factor LT
Matrix Cholesky::get_LT()
{
  assert(isDecomposed == true);
  switch(Decompose) {
      case chol_STANDARD:
         return (*LT1);
         break;
      case chol_GENERAL:
         return (*LT2);
         break;
  }
  //--- if decomposition is NOT performed
  cerr << "\nDecomposition MUST be done first!\n";
  return A;
}

//! get diagonal elements of decomposition (for A=L*LT or A=L*D*LT)
Vector Cholesky::get_D()
{
  assert(isDecomposed == true);
  switch(Decompose) {
      case chol_STANDARD:
         return D1;
         break;
      case chol_GENERAL:
         return D2;
         break;
  }
  //--- if decomposition is NOT performed
  cerr << "\nDecomposition MUST be done first!\n";
  Vector dd(N);
  return dd;
}

//! get inverse L_inv of lower Cholesky factor L
Matrix Cholesky::get_L_inv()
{
  if (isInverse == true)
    return (*L_inv);
  //--- if decomposition and inversion are NOT performed
  cerr << "\nDecomposition MUST be done first!\n";
  cerr << "\nAfter that, call Inversion!\n";
  return A;
}

//! get determinant of the input matrix A
const double Cholesky::getDeterminant() const
{
  if (isDecomposed == true)
    return Det;
  //--- if decomposition is NOT performed
  cerr << "\nDecomposition MUST be done first!\n";
  return (double)1.0;
}

//! get inverse of the input matrix A
Matrix Cholesky::getInverse()
{
  if (isInverse == true)
    return (*A_inv);
  //--- if decomposition and inversion are NOT performed
  cerr << "\nDecomposition MUST be done first!\n";
  cerr << "\nAfter that, call Inversion!\n";
  return A;
}

//! get solution of linear equations A * x = b
Vector Cholesky::getSolution()
{
  assert(isSolved == true);
  return x;
}

//! get one solutions of the system of linear equations A * X = B
Vector Cholesky::getSolutionNo(int i)
{
  assert(isSolved == true);
  assert(i >=0 && i<M);
  return (X.getColData(i));
}

//! get all solutions of the system of linear equations A * X = B
Matrix Cholesky::getAllSolutions()
{
  assert(isSolved == true);
  return X;
}

//========================================== private working functions
//! Standard Cholesky decomposition of SPD matrix A: A = L * LT
void Cholesky::cholDec()
{
  int i,j,k;
  double sum, temp, zero = 0.0;
  for (i=0; i<N; i++) {
    for (j=i; j<N; j++) {
      sum = A.getElement(i,j);
      for (k=0; k<i; k++) {
        sum -= L1.getElement(i,k) * L1.getElement(j,k);
      }
      if (i == j) {
        if (sum <= zero) {
          // it should never be here
          cout << "\nCHOLDEC(): Matrix is NOT positive definite!" << endl;
          return;
        }
        D1.setData(i,sqrt(sum));   // diagonal elements of L are in vector D
      } // end if i == j
      else {
        temp = sum / D1.getData(i);
        L1.setElement(j,i, temp);  // lower part of the matrix A = L
      }
    } // end for j
  } // end for i

  //--- put diagonal elements into position
  for (i = 0; i < N; i++)
    L1.setElement(i,i, D1.getData(i));
  //--- give the name
  L1.setName("Lower_Triangle");
  //--- make the other cholesky factor
  LT1 = new Matrix(N,N);
  for (i=0; i<N; i++)
    for (j=i; j<N; j++)
      LT1->setElement(i,j,L1.getElement(j,i));
  LT1->setName("Upper_Triangle");
}

//! Generalized Cholesky decomposition of SPD matrix A: A = L * D * LT
void Cholesky::cholDecGen()
{
  int i,j,k;
  double sum, temp, zero = 0.0;
  for (i=0; i<N; i++) {
    for (j=i; j<N; j++) {
      sum = A.getElement(i,j);
      for (k=0; k<i; k++) {
        sum -= L2.getElement(i,k) * L2.getElement(j,k);
      }
      if (i == j) {
        if (sum <= zero) {
          cout << "\nCHOLDEC(): Matrix is NOT positive definite!" << endl;
        }
        D2.setData(i,sqrt(sum));   // diagonal elements of L are in vector D
      } // end if i == j
      else {
        temp = sum / D2.getData(i);
        L2.setElement(j,i, temp);  // lower part of the matrix A = L
      }
    } // end for j
  } // end for i
  double di, Lji;
  for (i=0; i<N-1; i++) {
    di = D2.getData(i);
    for (j=i+1; j<N; j++) {
      Lji = L2.getElement(j,i);
      L2.setElement(j,i, Lji/di);
    }
  }
  //--- put diagonal elements into position and form vector D
  for (i = 0; i < N; i++) {
    L2.setElement(i,i, 1.0);
    di = D2.getData(i);
    D2.setData(i,di*di);
  }
  //--- give the name
  L2.setName("Lower_Triangle");
  //--- make the other cholesky factor
  LT2 = new Matrix(N,N);
  for (i=0; i<N; i++)
    for (j=i; j<N; j++)
      LT2->setElement(i,j,L2.getElement(j,i));
  LT2->setName("Upper_Triangle");
}

//! inverse of a matrix A
Matrix Cholesky::cholInv()
{
  //--- Matrix L_inv - inverse of Cholesky factor L
  L_inv = new Matrix(N,N);
  assert(L_inv != NULL);
  int i,j,k;
  double sum;
  for (i=0; i<N; i++) {
    L_inv->setElement(i,i, 1.0/L1.getElement(i,i));
    for (j=i+1; j<N; j++) {
      sum = 0.0;
      for (k=i; k<j; k++)
        sum -= L1.getElement(j,k) * L_inv->getElement(k,i);
      L_inv->setElement(j,i, sum/D1.getData(j));
    }
  }
  //--- give the name
  L_inv->setName("Inverse_L");
  //--- obtain the inverse of the input matrix A
  double Li,Lj,Ai,Aij;
  Matrix ans(N,N);               // inverse matrix is ans

  for (i=0; i<N; i++) {
    Li = L_inv->getElement(i,i);
    ans.setElement(i,i, Li*Li);
    for (k=i+1; k<N; k++) {
      Li = L_inv->getElement(k,i);
      Ai = ans.getElement(i,i);
      Ai += Li*Li;
      ans.setElement(i,i, Ai);
    }
    for (j=i+1; j<N; j++) {
      for (k=j; k<N; k++) {
        Aij = ans.getElement(i,j);
        Li  = L_inv->getElement(k,i);
        Lj  = L_inv->getElement(k,j);
        Aij += Li*Lj;
        ans.setElement(i,j, Aij);
      }
    } // end for j=i+1
  } // end for i
  for (i=0; i<N; i++) {
    for (j=0; j<i; j++) {
      ans.setElement(i,j, ans.getElement(j,i));
    }
  }
  //--- give the name
  ans.setName("Inverse_A");
  isInverse = true;
  return (ans);
}


//! solution of linear equations A * x = b
Vector Cholesky::cholEqs(const Vector& b)
{
  //--- result vector x
  Vector ans(N);
  int i,k;
  double sum, L_ik, x_k, D_i;
  for (i=0; i<N; i++) {         // forward solution L * y = b
    D_i = D1.getData(i);
    sum = b.getData(i);
    for (k=i-1; k>=0; k--) {
      L_ik = L1.getElement(i,k);
      x_k  = ans.getData(k);
      sum -= L_ik * x_k;
    }
    ans.setData(i, sum/D_i);
  }
  for (i=N-1; i>=0; i--) {      // backsubstitution LT * x = y
    D_i = D1.getData(i);
    sum = ans.getData(i);
    for (k=i+1; k<N; k++) {
      L_ik = L1.getElement(k,i);
      x_k  = ans.getData(k);
      sum -= L_ik * x_k;
    }
    ans.setData(i, sum/D_i);
  }
  ans.setName("Solution");
  return (ans);
}

//! solution of system of linear equations A * X = B
Matrix Cholesky::cholEqs(const Matrix& BB)
{
  //--- result matrix
  Matrix ans(N,M);
  //--- right-hand side vector no. i
  Vector* bi = new Vector(N);
  assert(bi != NULL);
  //--- solution vector no. i
  Vector* xi = new Vector(N);
  assert(xi != NULL);

  string solName = "Solution No. ";
  string broj;
  char Broj[10];
  //--- now, solve equations M times
  for (int ii=0; ii<M; ii++) {
    //--- obtain solution
    *bi = BB.getColData(ii);
    *xi = cholEqs(*bi);
    sprintf(Broj, "%d", (ii+1));
    broj = Broj;
    solName += broj;
    xi->setName(solName);
    //--- transfer solution
    ans.setColData(ii,*xi);
  }
  // release dynamic memory
  delete xi;
  delete bi;
  // return solution
  return (ans);
}

//! control of Cholesky decomposition
Matrix Cholesky::cholControl()
{
  Matrix res(N,N);
  switch (Decompose) {
      case chol_STANDARD:
         res = L1*(*LT1);
         break;
      case chol_GENERAL:
         Matrix* Dia = new Matrix(N,N);
         assert(Dia != NULL);
         Dia->Diagonal(D2);
         res = L2*(*Dia)*(*LT2);
         delete Dia;
         break;
  }
  return (res);
}

}   // end of namespace smv_lib

//================================ end of file
